{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Effects of different loss functions for classification\n",
    "In this jupyter notebook assignment, we will investigate the use of different methods to do binary classification and see how things work out. In particular, we will investigate four different approaches for learning classification boundaries:\n",
    "1. Mean classification, i.e., calculating means of the labeled training data and classify based on distances from these means\n",
    "2. Ordinary least squares, i.e., we calculate a best linear line that separates data based on standard squared loss\n",
    "3. Exponential regression, i.e., we calculate a best linear line that separates data based on exponential loss\n",
    "4. logistic regression, i.e., we calculate a best linear line that separates data based on logistic loss --- a \"mellower\" version of exponential loss that is less strict in penalizing \"big\" classification errors.\n",
    "\n",
    "More details will be included while we go through the examples. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first thing that we will do is simply illustrate the losses 1,2,3 above so you can see how they look different. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_x = np.arange(-2, 2, 0.01)\n",
    "squared_error = [(i-1)**2 for i in line_x]\n",
    "squared_error_minus = [(i+1)**2 for i in line_x]\n",
    "exp_error = np.exp(line_x)\n",
    "exp_error_minus = np.exp(-line_x)\n",
    "logistic_error = np.log(1+np.exp(line_x))\n",
    "logistic_error_minus = np.log(1+np.exp(-line_x))\n",
    "abs_val = np.abs(line_x)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,10))\n",
    "ax.grid(True, which='both')\n",
    "plt.plot(line_x, squared_error, '-.', label='Squared Error for - points', color='r')\n",
    "plt.plot(line_x, exp_error, '-.', label='Exponential Error for - points', color='b')\n",
    "plt.plot(line_x, logistic_error, '-.', label='Logistic Error for - points', color='g')\n",
    "plt.plot(line_x, squared_error_minus, label='Squared Error for + points', color='r')\n",
    "plt.plot(line_x, exp_error_minus, label='Exponential Error for + points', color='b')\n",
    "plt.plot(line_x, logistic_error_minus, label='Logistic Error for + points', color='g')\n",
    "plt.plot(line_x, abs_val, ':', label='y = |x|', color='k')\n",
    "plt.title('Three Different Errors')\n",
    "plt.legend()\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "ax.set_ylim([-0.25,4])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this plot, we plot out the error shapes for the $+1$ and $-1$ labeled classes. Notice that the logistic error approaches the plot of $y = |x|$ on one side and the exponential error on the other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating dataset\n",
    "\n",
    "We generate a synthetic dataset by selecting red and blue points such that they lie approximately in a line.\n",
    "\n",
    "For our different examples, we will be generating the training data in different ways to illustrate the capability of each approach in dealing with challenges. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_generate_points(parameters, num_points, outlier_points, extreme_outlier_points, red = True, example4 = False):\n",
    "    \"\"\"\n",
    "    Randomly generate points such that a large cluster lives at [1, a+b] for num_points-outlier_points points, and the rest lives on the line ax+b with noise.\n",
    "    args:\n",
    "        parameters: a length 3 list [a b sigma] that describes the underlying true linear line y = ax+b and the level of loss\n",
    "        num_points: number of clustered data points near the true line\n",
    "        outlier_points: number of non-clustered datapoints near the true line\n",
    "        extreme_outlier_points: number of datapoints far from true line\n",
    "        red: a boolean variable signaling if generated data points are red points (as opposed to blue points)\n",
    "        example4: a boolean variable signaling if we are generating points for Example 4. \n",
    "    \"\"\"\n",
    "    a = parameters[0]\n",
    "    b = parameters[1]\n",
    "    sigma = parameters[2]\n",
    "    \n",
    "    x1 = np.random.uniform(0.8,1.2, size = [num_points,1])\n",
    "    x2 = np.random.uniform(-5, 5, size = [outlier_points,1])\n",
    "    x = np.concatenate((x1, x2))\n",
    "    points = np.random.normal(a*x+b, sigma, size = [num_points+ outlier_points,1])\n",
    "    out = np.concatenate((x,points),axis = 1)   \n",
    "    x3 = np.random.uniform(2, 5, size = [extreme_outlier_points,1])\n",
    "    \n",
    "    if not example4:\n",
    "        if red == True:\n",
    "            extreme_points = np.random.normal(a*x3+b + 20, sigma, size = [extreme_outlier_points,1])\n",
    "        else:\n",
    "            x3 = np.random.uniform(-30, -25, size = [extreme_outlier_points,1])\n",
    "            extreme_points = np.random.normal(a*x3+b - 50, sigma, size = [extreme_outlier_points,1])\n",
    "        ext = np.concatenate((x3, extreme_points),axis = 1)\n",
    "\n",
    "        out = np.concatenate((out, ext))\n",
    "\n",
    "        return out\n",
    "    else:\n",
    "        x4 = np.random.uniform(10, 12, size = [extreme_outlier_points,1])\n",
    "        if red == True:\n",
    "            extreme_points = np.random.normal(a*x4+b, sigma, size = [extreme_outlier_points,1])\n",
    "        else:\n",
    "            extreme_points = np.random.normal(a*x4+b+30, 10*sigma, size = [extreme_outlier_points,1])\n",
    "        ext = np.concatenate((x4, extreme_points),axis = 1)\n",
    "\n",
    "        out = np.concatenate((out, ext))\n",
    "\n",
    "        return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1: Mean Classification\n",
    "\n",
    "In this example, we showcase a simple example where mean classification performs nicely. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_points_parameters = [2, 5, 1] # = [a,b,sigma], assuming y ~ ax+b \n",
    "blue_points_parameters = [-2, -5, 1]\n",
    "\n",
    "num_points = 100\n",
    "\n",
    "ar_red = random_generate_points(red_points_parameters, num_points, 0, 0)\n",
    "ar_blue = random_generate_points(blue_points_parameters, num_points, 0, 0, red = False)\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Training Data')\n",
    "plt.legend()\n",
    "ax.set_ylim([-10,10])\n",
    "ax.set_xlim([0,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this set of training data, we can calculate the means of the red and blue points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_c = np.mean(ar_red, axis=0)\n",
    "blue_c = np.mean(ar_blue, axis=0)\n",
    "coeffs = red_c - blue_c\n",
    "mid = (red_c + blue_c)/2.0\n",
    "\n",
    "line_x = np.linspace(-5, 5)\n",
    "line_y2 = [-(coeffs[0]/coeffs[1])*(i - mid[0]) + mid[1] for i in line_x]\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=3, alpha = 0.3)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=3, alpha = 0.3)\n",
    "plt.scatter(red_c[0], red_c[1], s=140, label='Red Center', marker='*', color='r')# facecolors='none', edgecolors='r')\n",
    "plt.scatter(blue_c[0], blue_c[1], s=140, label='Blue Center', marker='*', color='b')# facecolors='none', edgecolors='r')\n",
    "plt.plot(line_x, line_y2, label='Classifier Boundary', color='k')\n",
    "plt.legend()\n",
    "ax.set_ylim([-10,10])\n",
    "ax.set_xlim([0,2])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2: Mean classification fails but ordinary least squares work\n",
    "\n",
    "In the above example, we saw that mean classification works nicely when the underlying distribution and the data points are very nice behaving. In this example, we will see a case when mean classification fails but ordinary least squares work.\n",
    "\n",
    "Again, we start by generating our training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_points_parameters = [-2, 5, 1] # = [a,b,sigma], assuming y ~ ax+b \n",
    "blue_points_parameters = [-2, -5, 1]\n",
    "\n",
    "num_points = 100 # number of red or blue points\n",
    "outlier_points = 10\n",
    "extreme_outlier_points = 20\n",
    "\n",
    "ar_red = random_generate_points(red_points_parameters, num_points, outlier_points, 0)\n",
    "ar_blue = random_generate_points(blue_points_parameters, num_points, outlier_points, 0, red = False)\n",
    "\n",
    "# Plot dataset\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Training points')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we calculate the ordinary least squares classification boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.concatenate((ar_red, ar_blue))\n",
    "n = len(X)\n",
    "\n",
    "# Add column of ones to account for bias term\n",
    "X = np.concatenate((X, np.ones((n,1))), axis=1)\n",
    "\n",
    "# Create y array of class labels\n",
    "y = np.concatenate((np.ones(int(n/2)), -np.ones(int(n/2))))\n",
    "\n",
    "# Calculate the Ordinary Least Squares solution\n",
    "beta = np.linalg.inv(X.T @ X) @ (X.T @ y)\n",
    "\n",
    "# Plot the resulting regression line\n",
    "line_x = np.linspace(-5, 5)\n",
    "# line_y = -beta[0] / beta[2] - (beta[1] / beta[2]) * line_x\n",
    "line_y = -beta[0] / beta[1] * line_x - beta[2] / beta[1] \n",
    "\n",
    "# Plot dataset + centroid + classification\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Ordinary Least Squares Classification Boundary')\n",
    "plt.plot(line_x, line_y, label='Classifier Boundary', color='k')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "red_c = np.mean(ar_red, axis=0)\n",
    "blue_c = np.mean(ar_blue, axis=0)\n",
    "coeffs = red_c - blue_c\n",
    "mid = (red_c + blue_c)/2.0\n",
    "line_y2 = [-(coeffs[0]/coeffs[1])*(i - mid[0]) + mid[1] for i in line_x]\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=1, alpha = 0.5)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=1, alpha = 0.5)\n",
    "plt.scatter(red_c[0], red_c[1], s=140, label='Red Center', marker='*', color='r')# facecolors='none', edgecolors='r')\n",
    "plt.scatter(blue_c[0], blue_c[1], s=140, label='Blue Center', marker='*', color='b')# facecolors='none', edgecolors='r')\n",
    "plt.title('Mean classification Boundary')\n",
    "plt.plot(line_x, line_y2, label='Classifier Boundary', color='k')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, ordinary least squares works very well for this classification task but mean-based classification is making some mistakes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 3: Exponential regression\n",
    "\n",
    "In this example, we will see an example where using exponential regression can be a good choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_points_parameters = [-2, 5, 2] # = [a,b,sigma], assuming y ~ ax+b \n",
    "blue_points_parameters = [-2, -5, 2]\n",
    "\n",
    "num_points = 10 # number of red or blue points\n",
    "outlier_points = 30\n",
    "extreme_outlier_points = 40\n",
    "\n",
    "ar_red = random_generate_points(red_points_parameters, num_points, outlier_points, extreme_outlier_points)\n",
    "ar_blue = random_generate_points(blue_points_parameters, num_points, outlier_points, extreme_outlier_points, red = False)\n",
    "\n",
    "# Plot dataset\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Training points')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we calculate the ordinary least squares classification boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.concatenate((ar_red, ar_blue))\n",
    "n = len(X)\n",
    "\n",
    "# Add column of ones to account for bias term\n",
    "X = np.concatenate((X, np.ones((n,1))), axis=1)\n",
    "\n",
    "# Create y array of class labels\n",
    "y = np.concatenate((np.ones(int(n/2)), -np.ones(int(n/2))))\n",
    "\n",
    "# Calculate the Ordinary Least Squares solution\n",
    "beta = np.linalg.inv(X.T @ X) @ (X.T @ y)\n",
    "\n",
    "# Plot the resulting regression line\n",
    "line_x = np.linspace(-30, 6)\n",
    "# line_y = -beta[0] / beta[2] - (beta[1] / beta[2]) * line_x\n",
    "line_y = -beta[0] / beta[1] * line_x - beta[2] / beta[1] \n",
    "\n",
    "# Plot dataset + classification\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Classification Boundary of Ordinary Least Squares')\n",
    "plt.plot(line_x, line_y, label='Classifier Boundary', color='k')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's calculate our exponential regression boundary. To do this, we will be handcrafting an exponential regression where we use an exponential loss --- we want + points (red ones) to be positive and - points (blue ones) to be negative. More positive is better for + points and more negative is better for - points. This is done by applying $e^{-p}$ as the cost for points that should have a positive label and $e^{+p}$ as the cost for points that should have a negative label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_weight_exponential(w, X0, X1):\n",
    "    \"\"\"\n",
    "    Calculates the next weight for exponential regression\n",
    "    args:\n",
    "    w: current weight\n",
    "    X0: red training data\n",
    "    X1: blue training data\n",
    "    \"\"\"\n",
    "    p0 = np.dot(X0,w)\n",
    "    p1 = np.dot(X1,w)\n",
    "    \n",
    "    # -------------------------------------------------------------\n",
    "    # In this TODO section, you need to calculate the derivatives and then update weight w. \n",
    "    # As an example, if we have a loss function as a sum of (p0[i]-1)^2 corresponding to the \n",
    "    # square loss for red points,\n",
    "    # The first derivative first_der0 would be a vector with each element being 2(p0[i]-1)\n",
    "    # The second dsecond_der0 would be a vector with each element being 2\n",
    "    # and from derivations in lecture and homework,\n",
    "    # new_targets_0 would be a vector with each element being -2(p0[i]-1)/sqrt(2)\n",
    "    # new_X0 would be a matrix with each element of X multiplied by sqrt(2)\n",
    "    # One way to implement this would be\n",
    "    #     first_der0 = 2*p0 - 2    \n",
    "    #     second_der0 = 2*np.ones(p0.shape)\n",
    "    #     new_targets_0 = -first_der0/np.sqrt(2)\n",
    "    #     new_X0 = X0 * np.sqrt(2)\n",
    "    # You can plug this in and see that this will give the same solution as least squares.\n",
    "    # -------------------------------------------------------------\n",
    "    \n",
    "    first_der0 = # TODO: Calculate first derivative of loss function of p0\n",
    "    second_der0 = # TODO: Calculate second derivative of loss function of p0\n",
    "    \n",
    "    first_der1 = # TODO: Calculate first derivative of loss function of p1\n",
    "    second_der1 = # TODO: Calculate second derivative of loss function of p1\n",
    "    \n",
    "    # Compute the new targets based on the first and second derivatives\n",
    "    new_targets_0 = # TODO: Calculate updated target y corresponding to the red training data\n",
    "    new_targets_1 = # TODO: Calculate updated target y corresponding to the blue training data\n",
    "    \n",
    "    # reweight the X features to reflect what the second derivative tells us\n",
    "    new_X0 = # TODO: Calculate updated X corresponding to the red training data\n",
    "    new_X1 = # TODO: Calculate updated X corresponding to the blue training data\n",
    "    \n",
    "    #--------------------------------------------------------------\n",
    "    \n",
    "    new_X = np.concatenate((new_X0, new_X1))\n",
    "   \n",
    "    # Join the targets together into one long vector.\n",
    "    new_y = np.concatenate((new_targets_0, new_targets_1))\n",
    "    \n",
    "    # Now return the updated weight --- which is the old weight plus the least-squares solution for the new problem\n",
    "    return w + np.linalg.inv(new_X.T @ new_X) @ (new_X.T @ new_y)\n",
    "    \n",
    "    \n",
    "w = np.zeros((3,1))\n",
    "concatenated_red = np.concatenate((ar_red, np.ones((ar_red.shape[0],1))), axis = 1)\n",
    "concatenated_blue = np.concatenate((ar_blue, np.ones((ar_blue.shape[0],1))), axis = 1)\n",
    "\n",
    "for i in range(1000):\n",
    "    w = update_weight_exponential(w, concatenated_red, concatenated_blue)\n",
    "    \n",
    "def get_line_val(x, w):\n",
    "    return (-w[0]/w[1]) * x - w[2] / w[1]\n",
    "    \n",
    "line_x_axis = np.arange(-30,6,1)\n",
    "line_y_axis = [get_line_val(i, w) for i in line_x_axis]\n",
    "\n",
    "true_y_axis = [red_points_parameters[0]*i + (red_points_parameters[1]+blue_points_parameters[1])/2 for i in line_x_axis]\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=3)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=3)\n",
    "plt.title('Classification Boundaries')\n",
    "plt.plot(line_x_axis, line_y_axis, label='Exp. Classifier Boundary', color='m')\n",
    "plt.plot(line_x_axis, true_y_axis, label='True Classifier Boundary', color='g')\n",
    "plt.plot(line_x, line_y, label='Lst Sq Classifier Boundary', color='k')\n",
    "ax.set_ylim([-20,25])\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 4: Logistic Regression\n",
    "\n",
    "In this example, we will handcraft a logistic regression and see examples where logistic regression performs better than exponential regression. Again, we start by generating training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_points_parameters = [-2, 8, 1] # = [a,b,sigma], assuming y ~ ax+b \n",
    "blue_points_parameters = [-2, 3, 1]\n",
    "\n",
    "num_points = 20 # number of red or blue points\n",
    "outlier_points = 30\n",
    "extreme_outlier_points = 1\n",
    "\n",
    "ar_red = random_generate_points(red_points_parameters, num_points, outlier_points, \n",
    "                                extreme_outlier_points, example4 = True)\n",
    "ar_blue = random_generate_points(blue_points_parameters, num_points, outlier_points, \n",
    "                                 extreme_outlier_points, red = False, example4 = True)\n",
    "\n",
    "# Plot dataset\n",
    "ax = plt.subplot(111)\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Training points')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get our weight estimations by using our weight updates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def update_weight_logistic(w, X0, X1):\n",
    "    \"\"\"\n",
    "    Calculates the next weight for exponential regression\n",
    "    args:\n",
    "    w: current weight\n",
    "    X0: red training data\n",
    "    X1: blue training data\n",
    "    \"\"\"\n",
    "    p0 = np.dot(X0,w)\n",
    "    p1 = np.dot(X1,w)\n",
    "    \n",
    "    # -------------------------------------------------------------\n",
    "    \n",
    "    first_der0 = # TODO: Calculate first derivative of loss function of p0\n",
    "    second_der0 = # TODO: Calculate second derivative of loss function of p0\n",
    "    \n",
    "    first_der1 = # TODO: Calculate first derivative of loss function of p1\n",
    "    second_der1 = # TODO: Calculate second derivative of loss function of p1\n",
    "    \n",
    "    # Compute the new targets based on the first and second derivatives\n",
    "    new_targets_0 = # TODO: Calculate updated target y corresponding to the red training data\n",
    "    new_targets_1 = # TODO: Calculate updated target y corresponding to the blue training data\n",
    "    \n",
    "    # reweight the X features to reflect what the second derivative tells us\n",
    "    new_X0 = # TODO: Calculate updated X corresponding to the red training data\n",
    "    new_X1 = # TODO: Calculate updated X corresponding to the blue training data\n",
    "    \n",
    "    #--------------------------------------------------------------\n",
    "    \n",
    "    new_X = np.concatenate((new_X0, new_X1))\n",
    "    new_y = np.concatenate((new_targets_0, new_targets_1))\n",
    "    \n",
    "    return w + np.linalg.inv(new_X.T @ new_X) @ (new_X.T @ new_y)\n",
    "    \n",
    "    \n",
    "w_log = np.zeros((3,1))\n",
    "w_exp = np.zeros((3,1))\n",
    "concatenated_red = np.concatenate((ar_red, np.ones((ar_red.shape[0],1))), axis = 1)\n",
    "concatenated_blue = np.concatenate((ar_blue, np.ones((ar_blue.shape[0],1))), axis = 1)\n",
    "\n",
    "for i in range(100):\n",
    "    w_log = update_weight_logistic(w_log, concatenated_red, concatenated_blue)\n",
    "    w_exp = update_weight_exponential(w_exp, concatenated_red, concatenated_blue)\n",
    "\n",
    "line_x_axis = np.arange(-6,12,1)\n",
    "line_y_axis_log = [get_line_val(i, w_log) for i in line_x_axis]\n",
    "line_y_axis_exp = [get_line_val(i, w_exp) for i in line_x_axis]\n",
    "true_y_axis = [red_points_parameters[0]*i + (red_points_parameters[1]+blue_points_parameters[1])/2 for i in line_x_axis]\n",
    "\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title('Classification Boundaries')\n",
    "plt.plot(line_x_axis, line_y_axis_log, label='Log. Classifier Boundary', color='m')\n",
    "plt.plot(line_x_axis, line_y_axis_exp, label='Exp. Classifier Boundary', color='k')\n",
    "plt.plot(line_x_axis, true_y_axis, label = 'True Classifier Boundary', color = 'g')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, even with one mislabeled outlier we are able to shift the exponential classification boundary by a large amount. Meanwhile, logistic regression is indeed more mellow about this. This illustrates how logistic regression is not as sensitive to outliers as the exponential regression model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's see how the classifier boundary changes for both exponential and logistic regression in each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_log = np.zeros((3,1))\n",
    "w_exp = np.zeros((3,1))\n",
    "\n",
    "X = np.concatenate((ar_red, ar_blue))\n",
    "n = len(X)\n",
    "\n",
    "# Add column of ones to account for bias term\n",
    "X = np.concatenate((X, np.ones((n,1))), axis=1)\n",
    "\n",
    "# Create y array of class labels\n",
    "y = np.concatenate((np.ones(int(n/2)), -np.ones(int(n/2))))\n",
    "\n",
    "# Calculate the Ordinary Least Squares solution\n",
    "beta = np.linalg.inv(X.T @ X) @ (X.T @ y)\n",
    "\n",
    "# Plot the resulting regression line\n",
    "line_y_LS = -beta[0] / beta[1] * line_x_axis - beta[2] / beta[1] \n",
    "\n",
    "\n",
    "for i in range(4):\n",
    "    w_log = update_weight_logistic(w_log, concatenated_red, concatenated_blue)\n",
    "    w_exp = update_weight_exponential(w_exp, concatenated_red, concatenated_blue)\n",
    "    #print(w_log, w_exp)\n",
    "    line_x_axis = np.arange(-6,12,1)\n",
    "    line_y_axis_log = [get_line_val(i, w_log) for i in line_x_axis]\n",
    "    line_y_axis_exp = [get_line_val(i, w_exp) for i in line_x_axis]\n",
    "    true_y_axis = [red_points_parameters[0]*i + (red_points_parameters[1]+blue_points_parameters[1])/2 for i in line_x_axis]\n",
    "    plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "    plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "    plt.title(f\"Comparison at Iteration {i+1}\")\n",
    "    plt.plot(line_x_axis, line_y_axis_exp, label = 'Exp. Classifier Boundary', color = 'm')\n",
    "    plt.plot(line_x_axis, line_y_axis_log, label = 'Log. Classifier Boundary', color = 'c')\n",
    "    plt.plot(line_x_axis, line_y_LS, label = 'Lst Sq Classifier Boundary', color = 'k')\n",
    "    plt.plot(line_x_axis, true_y_axis, label = 'True Classifier Boundary', color = 'g')\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that the line is almost fixed after a very small number of iterations. These iterated least-squares approaches are very fast to converge in general. This is why in practice, logistic or exponential regression costs almost the same to run as ordinary least squares."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing with sklearn's logistic regression function\n",
    "In this part we will see that our logistic regression gives the same solution as scikitlearn's logistic regression function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "clf = LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=1000000.0, \n",
    "                         fit_intercept=False, intercept_scaling=1, class_weight=None, \n",
    "                         random_state=None, solver='lbfgs', max_iter=100, \n",
    "                         multi_class='ovr', verbose=0, warm_start=False, n_jobs=None).fit(X, y)\n",
    "coeffs = clf.coef_[0]\n",
    "line_y_axis_sklearn = [get_line_val(i, coeffs) for i in line_x_axis]\n",
    "line_y_axis_log = [get_line_val(i, w_log) for i in line_x_axis]\n",
    "\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title(f\"Sklearn's logistic boundary\")\n",
    "plt.plot(line_x_axis, line_y_axis_sklearn, label = 'sklearn log. Classifier Boundary', color = 'y')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(ar_red[:,0], ar_red[:,1], label='red', color='r', s=10)\n",
    "plt.scatter(ar_blue[:,0], ar_blue[:,1], label='blue', color='b', s=10)\n",
    "plt.title(f\"Our logistic boundary\")\n",
    "plt.plot(line_x_axis, line_y_axis_log, label = 'Log. Classifier Boundary', color = 'c')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "print(f\"Function of sklearn boundary: {coeffs[0]}x + {coeffs[1]}y + {coeffs[2]} = 0.\")\n",
    "print(f\"Function of our boundary: {w_log[0][0]}x + {w_log[1][0]}y + {w_log[2][0]} = 0.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that our logistic regression boundary is practically the same as sklearn's boundary, with very small numerical difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}